PCA analysis of extracted metrics - updated to Noah’s full set of data
ACS
helper function to process time series
And some package and labeling info
library(ggplot2)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)
library(factoextra) # easier visualizing outputs of prcomp## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBatimeseries_dir <- 'extracted_timeseries/'
metrics_write_dir <- 'extracted_timeseries/extracted_metrics/'
# get ecs ordering/labels
esm_labels <- read.csv(paste0(timeseries_dir,'global_tas_allesms.csv'), stringsAsFactors = FALSE) %>%
select(esm) %>% distinct %>%
mutate(plotesm = paste0(letters[as.integer(row.names(.))], '.', esm),
ECS_order = as.integer(row.names(.)))
print(esm_labels)## esm plotesm ECS_order
## 1 CAMS-CSM1-0 a.CAMS-CSM1-0 1
## 2 MIROC6 b.MIROC6 2
## 3 GFDL-ESM4 c.GFDL-ESM4 3
## 4 FGOALS-g3 d.FGOALS-g3 4
## 5 MPI-ESM1-2-HR e.MPI-ESM1-2-HR 5
## 6 MPI-ESM1-2-LR f.MPI-ESM1-2-LR 6
## 7 MRI-ESM2-0 g.MRI-ESM2-0 7
## 8 ACCESS-ESM1-5 h.ACCESS-ESM1-5 8
## 9 IPSL-CM6A-LR i.IPSL-CM6A-LR 9
## 10 CESM2-WACCM j.CESM2-WACCM 10
## 11 UKESM1-0-LL k.UKESM1-0-LL 11
## 12 CanESM5 l.CanESM5 12Spatial info
we want the shapes for plotting anyway, so prep them
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE## Reading layer `IPCC-WGI-reference-regions-v4' from data source
## `/Users/snyd535/Documents/GitHub/stitches_in_r/R/inst/shinyApp/python_curation/IPCC-WGI-reference-regions-v4_shapefile/IPCC-WGI-reference-regions-v4.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 58 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 90
## Geodetic CRS: WGS 84
Load ESM data that’s been processed
region_summary_main <- read.csv(paste0(metrics_write_dir, 'IPCC_land_regions_metrics.csv'),
stringsAsFactors = FALSE)
region_summary_main %>%
filter(experiment != 'ssp119',
experiment != 'ssp434',
experiment != 'ssp460',
experiment != 'ssp534-over') %>%
rename(region = acronym) ->
region_summary
print(head(region_summary))## esm experiment ensemble variable type name region
## 1 ACCESS-CM2 ssp126 r1i1p1f1 pr Land Arabian-Peninsula ARP
## 2 ACCESS-CM2 ssp126 r1i1p1f1 pr Land Central-Africa CAF
## 3 ACCESS-CM2 ssp126 r1i1p1f1 pr Land-Ocean Caribbean CAR
## 4 ACCESS-CM2 ssp126 r1i1p1f1 pr Land C.Australia CAU
## 5 ACCESS-CM2 ssp126 r1i1p1f1 pr Land C.North-America CNA
## 6 ACCESS-CM2 ssp126 r1i1p1f1 pr Land E.Asia EAS
## iasd end_val end_anomaly end_anomaly_pct mid_century_val
## 1 7.227917e-07 1.700625e-06 4.306991e-07 0.339153029 2.051770e-06
## 2 3.009324e-06 4.828026e-05 -4.965491e-08 -0.001027416 4.941042e-05
## 3 6.106140e-06 3.052769e-05 2.089845e-07 0.006892922 3.429079e-05
## 4 2.830779e-06 1.130480e-05 -1.269868e-06 -0.100986213 1.126320e-05
## 5 2.861294e-06 3.232150e-05 1.953880e-06 0.064340892 3.222543e-05
## 6 2.687448e-06 5.449186e-05 5.826705e-06 0.119730531 5.260652e-05
## mid_anomaly mid_anomaly_pct
## 1 7.818446e-07 0.61566180
## 2 1.080502e-06 0.02235680
## 3 3.972086e-06 0.13101110
## 4 -1.311461e-06 -0.10429388
## 5 1.857811e-06 0.06117738
## 6 3.941364e-06 0.08098944Prep data for PCA
idea is that there’s something really similar across ESMs, especially for SSP126/245 -> their RGB maps have so much in common.
is there a sort of core map that gets deviations added to as you move through scenarios and ESMs? Or how many distinct maps do we actually have
And the aim of the archive is to establish that our selection of ESMs forms a spanning set of vectors for all/most of CMIP6
observation = scenarioXesm
individual variable is a regionXmetric
so each observation vector is
3*Nregionslong:<Tr1...Trn, Pr1...Prn, IASDr1...IASDrn >for a scenarioXesm <-> there are 48 observations for each of 129 variablesAnd then becaues we’re interested in a pseudo-basis for the CMIP6 archive, we want the training data to be all scenarios and ESMs
Out of sample tests -> how does basis do on overshoots that it doesn’t train on; withhold one of the 12 ESMs and learn from other 11 to see how it goes; then see how a completely external ESM does on training from all 12
nothing a priori to preclude us from doing EOFs on the gridded time series data, or the gridded metrics, but having on IPCC regions avoids the issue of all ESMs being on their own grids
and overall doing region metrics helps keep data size operating on manageable
# reshape:
grouped_data <- split(region_summary, list(region_summary$esm, region_summary$experiment))
grouped_data <- grouped_data[sapply(grouped_data, function(x) dim(x)[1]) > 0]
shaped_data <- lapply(grouped_data, function(group){
if (nrow(group) >0) {
group %>%
filter(variable == 'tas') %>%
select(esm, experiment, ensemble,type, region,
iasd, end_anomaly, mid_anomaly) %>%
rename(tas_iasd = iasd,
tas_end_anomaly = end_anomaly,
tas_mid_anomaly = mid_anomaly) %>%
left_join(group %>%
filter(variable == 'pr') %>%
select(esm, experiment, ensemble,type, region,
iasd, end_anomaly_pct, mid_anomaly_pct) %>%
rename(pr_iasd = iasd,
pr_end_anomaly_pct = end_anomaly_pct,
pr_mid_anomaly_pct = mid_anomaly_pct),
by = c('esm','experiment', 'ensemble', 'type', 'region')
) ->
reshaped
reshaped %>%
group_by(esm, experiment, type, region) %>%
summarize(tas_iasd=mean(tas_iasd),
tas_end_anomaly=mean(tas_end_anomaly, na.rm = T),
tas_mid_anomaly = mean(tas_mid_anomaly, na.rm = T),
pr_iasd = mean(pr_iasd, na.rm = T),
pr_end_anomaly_pct = mean(pr_end_anomaly_pct, na.rm = T),
pr_mid_anomaly_pct = mean(pr_mid_anomaly_pct, na.rm = T) ) %>%
ungroup() %>%
mutate(ensemble = 'ensemble_avg') -> #%>%
#bind_rows(reshaped) ->
reshaped2
#TODO need to change shaping if including ensemble members
reshaped2 %>%
select(-type) %>%
gather(metric, value, -esm, -experiment, -ensemble, -region) %>%
mutate(row_id = paste0(region, '~', metric),
col_id = paste0(esm, '~', experiment, '~', ensemble)) %>%
select(-esm, -experiment,-ensemble, -region, -metric) %>%
as.data.frame() ->
reshaped3
colnames(reshaped3) <- c(paste0(unique(reshaped3$col_id)), 'row_id', 'col_id')
rownames(reshaped3) <- paste0(reshaped3$row_id)
reshaped3 %>%
select(-col_id, -row_id) ->
out
return(out)
}
}
)## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'esm', 'experiment', 'type'. You can
## override using the `.groups` argument.# combine columns but then transpose because prcomp wants rows to be observations
# and columns to be variables (but doing it the other way first is easier to code)
data <- t(do.call(cbind, shaped_data) )
print(str(data))## num [1:111, 1:258] 0.331 0.291 0.278 0.323 0.311 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:111] "ACCESS-CM2~ssp126~ensemble_avg" "ACCESS-ESM1-5~ssp126~ensemble_avg" "BCC-CSM2-MR~ssp126~ensemble_avg" "CAMS-CSM1-0~ssp126~ensemble_avg" ...
## ..$ : chr [1:258] "ARP~tas_iasd" "CAF~tas_iasd" "CAU~tas_iasd" "CNA~tas_iasd" ...
## NULLdo PCA
have to remove some NA’s:
- manually checked all 4 - no PR ens average values
data1 <- na.omit(data)
print('removed rows due to missing data')## [1] "removed rows due to missing data"print(row.names(data)[c(which(!(row.names(data) %in% row.names(data1))))])## [1] "CIESM~ssp126~ensemble_avg" "CIESM~ssp245~ensemble_avg"
## [3] "CIESM~ssp585~ensemble_avg" "NorESM2-LM~ssp585~ensemble_avg"data_pca <- prcomp(data1, center=TRUE, scale = TRUE)
# quick visualize
factoextra::fviz_eig(data_pca)
further visualization
# coordinates of each ESM-SSP combo in the PCs
data_pca$x %>%
as.data.frame() %>%
mutate(id = row.names(.)) %>%
separate(id, into=c('model', 'scenario', 'ensemble'), sep = '~') %>%
select(model, scenario, PC1, PC2, PC3, PC4, PC5, PC6, PC7, PC8, PC9)->
coordinatesPC1 vs PC2
ggplot(coordinates) +
geom_point(mapping = aes(x = PC1, y = PC2, color = model, shape = scenario)) 
ggplot(coordinates) +
geom_text(mapping = aes(x = PC1, y = PC2, color = model, label = model), size=0.95) 
ggplot(coordinates) +
geom_point(mapping = aes(x = PC1, y = PC2, color = scenario))
PC1 vs PC3
ggplot(coordinates) +
geom_point(mapping = aes(x = PC1, y = PC3, color = model, shape = scenario))
ggplot(coordinates) +
geom_point(mapping = aes(x = PC1, y = PC3, color = scenario))
PC1 vs PC4
ggplot(coordinates) +
geom_point(mapping = aes(x = PC1, y = PC4, color = model, shape = scenario))
ggplot(coordinates) +
geom_point(mapping = aes(x = PC1, y = PC4, color = scenario))
PC2 vs PC3
ggplot(coordinates) +
geom_point(mapping = aes(x = PC2, y = PC3, color = model, shape = scenario))
ggplot(coordinates) +
geom_point(mapping = aes(x = PC2, y = PC3, color = scenario))
PC2 vs PC4
ggplot(coordinates) +
geom_point(mapping = aes(x = PC2, y = PC4, color = model, shape = scenario))
ggplot(coordinates) +
geom_point(mapping = aes(x = PC2, y = PC4, color = scenario))
PC3 vs PC4
ggplot(coordinates) +
geom_point(mapping = aes(x = PC3, y = PC4, color = model, shape = scenario))
ggplot(coordinates) +
geom_point(mapping = aes(x = PC3, y = PC4, color = scenario))
plot PCs as maps
pcs <- data_pca$rotation
str(pcs) ## num [1:258, 1:107] 0.0179 0.0129 0.0288 0.0218 0.0151 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:258] "ARP~tas_iasd" "CAF~tas_iasd" "CAU~tas_iasd" "CNA~tas_iasd" ...
## ..$ : chr [1:107] "PC1" "PC2" "PC3" "PC4" ...# pull the region and metric info back in
as.data.frame(pcs) %>%
mutate(row_id = row.names(.)) %>%
separate(row_id, into = c('region', 'metric'), sep ='~') ->
tmp
# We want to convert each set of 3 metrics in a region to a single
# rgb coded value.
# and we want all on same color scale.
# so reshape for each metric (which goes to one of r, g, b)
# with all PCs want to consider stacked.
tmp %>%
filter(metric == 'tas_end_anomaly') %>%
select(region, metric, PC1) %>%
rename(tas_end_anomaly = PC1) %>%
mutate(pc = 'PC1') %>%
bind_rows(tmp %>%
filter(metric == 'tas_end_anomaly') %>%
select(region, metric, PC2) %>%
rename(tas_end_anomaly = PC2) %>%
mutate(pc = 'PC2'),
tmp %>%
filter(metric == 'tas_end_anomaly') %>%
select(region, metric, PC3) %>%
rename(tas_end_anomaly = PC3) %>%
mutate(pc = 'PC3') ,
tmp %>%
filter(metric == 'tas_end_anomaly') %>%
select(region, metric, PC4) %>%
rename(tas_end_anomaly = PC4) %>%
mutate(pc = 'PC4'),
tmp %>%
filter(metric == 'tas_end_anomaly') %>%
select(region, metric, PC5) %>%
rename(tas_end_anomaly = PC5) %>%
mutate(pc = 'PC5'),
tmp %>%
filter(metric == 'tas_end_anomaly') %>%
select(region, metric, PC6) %>%
rename(tas_end_anomaly = PC6) %>%
mutate(pc = 'PC6'),
tmp %>%
filter(metric == 'tas_end_anomaly') %>%
select(region, metric, PC7) %>%
rename(tas_end_anomaly = PC7) %>%
mutate(pc = 'PC7'),
tmp %>%
filter(metric == 'tas_end_anomaly') %>%
select(region, metric, PC8) %>%
rename(tas_end_anomaly = PC8) %>%
mutate(pc = 'PC8'),
tmp %>%
filter(metric == 'tas_end_anomaly') %>%
select(region, metric, PC9) %>%
rename(tas_end_anomaly = PC9) %>%
mutate(pc = 'PC9'),
tmp %>%
filter(metric == 'tas_end_anomaly') %>%
select(region, metric, PC10) %>%
rename(tas_end_anomaly = PC10) %>%
mutate(pc = 'PC10')
) ->
tas
print(head(tas))## region metric tas_end_anomaly pc
## ARP~tas_end_anomaly...1 ARP tas_end_anomaly 0.09396300 PC1
## CAF~tas_end_anomaly...2 CAF tas_end_anomaly 0.09282450 PC1
## CAU~tas_end_anomaly...3 CAU tas_end_anomaly 0.09292208 PC1
## CNA~tas_end_anomaly...4 CNA tas_end_anomaly 0.09335706 PC1
## EAS~tas_end_anomaly...5 EAS tas_end_anomaly 0.09580033 PC1
## EAU~tas_end_anomaly...6 EAU tas_end_anomaly 0.09210200 PC1print(tail(tas))## region metric tas_end_anomaly pc
## WNA~tas_end_anomaly...425 WNA tas_end_anomaly 0.048192076 PC10
## WSAF~tas_end_anomaly...426 WSAF tas_end_anomaly 0.004906223 PC10
## WSB~tas_end_anomaly...427 WSB tas_end_anomaly 0.047285891 PC10
## CAR~tas_end_anomaly...428 CAR tas_end_anomaly 0.052966268 PC10
## MED~tas_end_anomaly...429 MED tas_end_anomaly 0.054653960 PC10
## SEA~tas_end_anomaly...430 SEA tas_end_anomaly 0.018720644 PC10tmp %>%
filter(metric == 'pr_end_anomaly_pct') %>%
select(region, metric, PC1) %>%
rename(pr_end_anomaly_pct = PC1) %>%
mutate(pc = 'PC1') %>%
bind_rows(tmp %>%
filter(metric == 'pr_end_anomaly_pct') %>%
select(region, metric, PC2) %>%
rename(pr_end_anomaly_pct = PC2) %>%
mutate(pc = 'PC2'),
tmp %>%
filter(metric == 'pr_end_anomaly_pct') %>%
select(region, metric, PC3) %>%
rename(pr_end_anomaly_pct = PC3) %>%
mutate(pc = 'PC3') ,
tmp %>%
filter(metric == 'pr_end_anomaly_pct') %>%
select(region, metric, PC4) %>%
rename(pr_end_anomaly_pct = PC4) %>%
mutate(pc = 'PC4') ,
tmp %>%
filter(metric == 'pr_end_anomaly_pct') %>%
select(region, metric, PC5) %>%
rename(pr_end_anomaly_pct = PC5) %>%
mutate(pc = 'PC5') ,
tmp %>%
filter(metric == 'pr_end_anomaly_pct') %>%
select(region, metric, PC6) %>%
rename(pr_end_anomaly_pct = PC6) %>%
mutate(pc = 'PC6') ,
tmp %>%
filter(metric == 'pr_end_anomaly_pct') %>%
select(region, metric, PC7) %>%
rename(pr_end_anomaly_pct = PC7) %>%
mutate(pc = 'PC7') ,
tmp %>%
filter(metric == 'pr_end_anomaly_pct') %>%
select(region, metric, PC8) %>%
rename(pr_end_anomaly_pct = PC8) %>%
mutate(pc = 'PC8') ,
tmp %>%
filter(metric == 'pr_end_anomaly_pct') %>%
select(region, metric, PC9) %>%
rename(pr_end_anomaly_pct = PC9) %>%
mutate(pc = 'PC9') ,
tmp %>%
filter(metric == 'pr_end_anomaly_pct') %>%
select(region, metric, PC10) %>%
rename(pr_end_anomaly_pct = PC10) %>%
mutate(pc = 'PC10')
) ->
pr
tmp %>%
filter(metric == 'tas_iasd') %>%
select(region, metric, PC1) %>%
rename(tas_iasd = PC1) %>%
mutate(pc = 'PC1') %>%
bind_rows(tmp %>%
filter(metric == 'tas_iasd') %>%
select(region, metric, PC2) %>%
rename(tas_iasd = PC2) %>%
mutate(pc = 'PC2'),
tmp %>%
filter(metric == 'tas_iasd') %>%
select(region, metric, PC3) %>%
rename(tas_iasd = PC3) %>%
mutate(pc = 'PC3') ,
tmp %>%
filter(metric == 'tas_iasd') %>%
select(region, metric, PC4) %>%
rename(tas_iasd = PC4) %>%
mutate(pc = 'PC4') ,
tmp %>%
filter(metric == 'tas_iasd') %>%
select(region, metric, PC5) %>%
rename(tas_iasd = PC5) %>%
mutate(pc = 'PC5'),
tmp %>%
filter(metric == 'tas_iasd') %>%
select(region, metric, PC6) %>%
rename(tas_iasd = PC6) %>%
mutate(pc = 'PC6'),
tmp %>%
filter(metric == 'tas_iasd') %>%
select(region, metric, PC6) %>%
rename(tas_iasd = PC6) %>%
mutate(pc = 'PC6'),
tmp %>%
filter(metric == 'tas_iasd') %>%
select(region, metric, PC7) %>%
rename(tas_iasd = PC7) %>%
mutate(pc = 'PC7'),
tmp %>%
filter(metric == 'tas_iasd') %>%
select(region, metric, PC8) %>%
rename(tas_iasd = PC8) %>%
mutate(pc = 'PC8'),
tmp %>%
filter(metric == 'tas_iasd') %>%
select(region, metric, PC9) %>%
rename(tas_iasd = PC9) %>%
mutate(pc = 'PC9'),
tmp %>%
filter(metric == 'tas_iasd') %>%
select(region, metric, PC10) %>%
rename(tas_iasd = PC10) %>%
mutate(pc = 'PC10')
) ->
tas_iasd
tas %>%
select(-metric) %>%
left_join(pr %>% select(-metric), by =c('region', 'pc')) %>%
left_join(tas_iasd %>% select(-metric), by =c('region', 'pc')) ->
colored_pcs
colored_pcs$r <- scales::rescale(colored_pcs$tas_end_anomaly, to =c(0,255))
colored_pcs$g <- scales::rescale(colored_pcs$tas_iasd, to =c(0,255))
colored_pcs$b <- scales::rescale(colored_pcs$pr_end_anomaly_pct, to =c(0,255))
colored_pcs %>%
left_join(as.data.frame(shp) %>% select(Acronym, mean_lon, mean_lat), by = c('region' = 'Acronym')) %>%
rename(lon = mean_lon, lat = mean_lat) %>%
mutate(color = rgb(.$r, .$g, .$b, maxColorValue = 255),
color_order = as.integer(row.names(.))) ->
pc_plotting
pc_plotting %>%
select(color_order, color) %>%
distinct() ->
colors
shp %>%
left_join(pc_plotting, by = c('Acronym' = 'region'))->
shp_pcs
p<- ggplot() +
geom_sf(data = shp_pcs %>% filter(pc != 'PC10') %>% na.omit, aes(fill = as.factor(color_order)) ) +
scale_fill_manual(values =colors$color)+
facet_wrap(~pc, ncol=3) +
theme(legend.position = 'none')
p
#look at whether Precip inc or dec -> not just the PCs' precip entries >= 0
# because standardized/
as.data.frame(data_pca$center) %>%
mutate(row_id = row.names(.)) %>%
filter(grepl('pr_end_anomaly_pct', row_id)) %>%
separate(row_id, into = c('region', 'trash'), sep = '~') %>%
select(-trash) ->
pr_centers
as.data.frame(data_pca$scale) %>%
mutate(row_id = row.names(.)) %>%
filter(grepl('pr_end_anomaly_pct', row_id)) %>%
separate(row_id, into = c('region', 'trash'), sep = '~') %>%
select(-trash) ->
pr_scales
shp_pcs %>%
left_join(pr_centers, by = c('Acronym'='region')) %>%
left_join(pr_scales, by = c('Acronym'='region')) %>%
mutate(precip_change = if_else( ((pr_end_anomaly_pct*`data_pca$scale`)- `data_pca$center`) >=0, 'inc', 'dec'))->
shp_pcs2
ggplot() +
geom_sf(data = shp_pcs2 %>% filter(pc != 'PC10') %>% na.omit, aes(#fill = as.factor(color_order),
color = precip_change )) +
#scale_fill_manual(values =colors$orig_color)+
scale_color_manual(values = c('red', 'blue'))+
facet_wrap(~pc, ncol=3) +
theme(legend.position = 'none') 